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Abstract 

A straightforward 3-point quadrature formula of closed type is de- 
rived that improves on Simpson's rule. Just using the additional in- 
formation of the integrand's derivative at the two endpoints we show 
the error is sixth order in grid spacing. Various error bounds for 
the quadrature formula are obtained to quantify more precisely the 
errors. Applications in numerical integration are given. With these 
error bounds, which are generally better than the usual Peano bounds, 
the composite formulas can be applied to integrands with lower order 
derivatives. 
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1 Introduction 

In recent years some authors have considered so called perturbed (corrected) 
quadrature rules. For example, the corrected midpoint and trapezoid quadra- 
ture rules are considered in [1] and [2]. As a specific introductory example, 
consider a correction to the midpoint rule. The classical midpoint rule has 
the form ^ 

f( x )dx = (b-a)f(^) + R(f) J (1) 



where R(f) is the remainder term (error) of this rule. A perturbed (corrected) 
rule involves the derivative at the endpoints: 

f f(x) dx=(b- a)f( a -f) + [/'(&) - f'(a)} + R,(f) , (2) 

where R\(f) is the remainder term (error) of this rule. The following prop- 
erties have analogues in the work we report on Simpson's rule: 

1. the original rule (1) is exact for polynomials of degree 1, while the 
perturbed rule is exact for polynomials of degree 3, 

2. estimates of the errors are respectively 

\R(f)\ < it ^ M ^ ( 3 ) 

(b - a) 3 

£ iwr M - (4) 

where \f"(t)\ < M 2 , t e [a, b] — the bound (4) is better than (3); 

3. a corresponding composite quadrature formula, for the corrected rule, 
has only one additional term, with respect to a composite formula for 
the original rule, 

4. the corrected composite formula has a better estimation of error than 
the original composite formula (a consequence of 2). 
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Indeed another bound is 

WW <!£(»-.)•, 

where \f (i) (t)\ < M 4 for t G [a, b] . 

The above properties are valid for similar corrected rules (for example, for 
the corrected trapezoid rule) [2] . However, we cannot correct all quadrature 
rules such that all the properties 1-4 hold. In Section 2 we show that the 
well-known Simpson's rule does have a simple endpoint correction, but that 
the quadrature weights have to be modified as well, see (8). We highlight 
some advantages of the corrected rule over the Simpson's rule. In Section 3 
various error bounds of this rule are obtained. These error bounds are gen- 
erally, but not always, better than the usual Peano error bounds. In Sec- 
tion 4 applications in numerical integration are given. An illustrative example 
demonstrates that the modified rule gives better results than Simpson's rule. 



2 Finite differences derive the modified 
Simpson's rule 

We modify Simpson's rule for integration. First, we consider integration 
rules formed over just two consectutive subintervals, each of length h, and 
derive endpoint modifications. This is analogous to the improvement to (1) 
achieved by the inclusion of endpoint derivative information in (2). Second, 
this modified Simpson's rule is straightforwardly summed to apply to an in- 
tegration over many subintervals. In later sections we rederive these formula 
with less restrictions on the integrand and with error bounds rather than just 
leading order estimates. 



Theorem 1 For C 1ac '[a, b] integrands f, 

ff(x)dx = b -^[7f( a ) + 16f(^) + 7f(b)] 
{b-af 



60 

where the error term is, to a leading order estimate 



[/'(&) - /'(a)] + R(f) , (5) 
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Example: As a simple illustrative example, consider f* e x dx = e — l/e = 
2.3504. Simpson's rule estimates the integral as approximately (e + 4 + 
l/e)/3 = 2.3621 , in error by about 0.01 . However, the modified Simpson's 
rule (5) estimates the integral as (6e + 16 + 8/e)/15 = 2.3502 , which has an 
error about two orders of magnitude smaller. This modification to Simpson's 
rule can be very effective. 

Proof: Consider integrating / over two consectutive intervals in a regular 
grid of points Xj with grid spacing h. Identify and hence 

Xj = (a + b)/2 is the midpoint. Following [3, p65], we write the analysis 
in terms of centred difference and mean operators, Sfj = fj+1/2 + fj-1/2 
and nfj = (/j+1/2 — f 3-1/2) /I respectively, and the differentiation operator 
denoted by d. Then the inter gral [3, p69] 



1 r x i+i 

— J f(x)dx = 0<9)~ V/j • 



To derive a three point integration rule with endpoint corrections such as (5), 
the above right-hand side must be in the form [1 + a5 2 + fi5/3(hd)]fj for some 
constant a and some function f3{hd): [1 + a5 2 ]fj symmetrically involves fj 
and fj±i alone; and /j,5/3(hd)fj only involves the derivatives of / at the end- 
points Xj±i . Thus we rearrange the operator equation 

(/id)" V =l+a5 2 + /i5/3(hd) 

no jj,o jj, 

& P(hd) = — - . I , _ - a2 tanh \hd , 
no sinn no 

as 5 = 2 sinh ^hd and fi = cosh ^hd [3, p65]. Choosing this particular func- 
tion /3(hd) would generate a rule in the interior together with end point 
corrections that would give an exact quadrature formula (for all a). How- 
ever, the infinite derivatives required are not practical. 

We chose a to generate an accurate rule that only needs to know func- 
tion values and the end point derivatives. The approach is to expand this 
function (3{hd) in powers of small hd to see 

P={-a+ \)hd + (i Q - ^){hdf + (-^a + ^){hdf + ■■■. 

Observe as an aside that choosing a = 1/6 eliminates the first derivative 
term in (3 leading to the familiar Simpson's rule with an error determined by 
the neglected parts of j3, namely the end-point contributions — fi5(hd) 3 / 180 + 
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■ • • - 1 Instead we choose a = 7/30 to eliminate the third derivative term in j3: 

I 3+1 f(x) dx = 2h(hd)- 1 / i5f j 

= 2M 1 + ^ fi + 2 M l-^hd + ^- Q hd 5 + •••]/, 
= ft [7/.^ + i 6/ . + 7fj+l] 

~ TE [fj+i ~~ f'j-i) + J25 [fj+i ~ fj-i] ^ • ( 7 ) 
Substitute /i = (6 — a) /2 to reproduce (5) and its leading order error (6). 4* 

Corollary 2 Apply (7) to n consecutive pairs of intervals from say a = xq 
to b = X2n and sum to immediately deduce that for C°°[a, b] integrands f 

r b h 2n-l ,2 

/ f(x)dx = - £ [7f J -i + lQf ] + 7f ]+1 ]--[f'(b)-f'(a)] + R(f). 

Ja j=l,j odd 

(8) 

where the error is to leading order 

h 6 

R(f)*^[f v (b)-na)} . (9) 

See that these simple modifications to Simpson's rule generate an inte- 
gration method with error (9(/i 6 ). 



3 Further error analysis 

For the sake of simplicity, we first consider the error analysis for the quadra- 
ture formula (5) on the interval [0,1]. Then we easily transform obtained 
results to an arbitrary interval. The formula (5) on the interval [0, 1] has the 
form 

// M ^ 7/(0)+16 3 f +^> -m_m +Mfh (10) 



lr rhus Simpson's rule has no simple end correction involving just the first derivative of 
the integrand, of the type described in (2) for the mid-point rule, instead the end correction 
would necessarily involve the third derivative. 
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for k = 2, 3, 4, 5, 6 . Using the Peano Kernel Theorem we find the following 
Peano kernels 

;-l) fc P fc (x), xe 
-l) k Q k (x), xe 



TJx) 



for k — 2, 3, 4, 5, 6 , where 



p 2 




— 1 ™2 7 ™ 1 1 
2"° 30 ^ 60 ' 


Q 2 




- 1^.2 _ 

2 


23 , 17 
30"° ~ 60 


p 3 




= |), 


Q 3 




= 31 (x — 


l) 2 (x-- 


Pa 




= ^x 2 (x-i)(x-|). 

— 51^ (% — 2)\ x ~ 3) • 


Q 4 




= h(?- 


P 5 




Q 5 




= |(*- 


l) 3 (x- - 
l) 4 (x 2 - 






_ 1 ™4/ 2 7 , 1\ 
gj-t' V A 5 ~ 2' ' 






— &\ x ~ 



25 > 
2^ 
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We have 



for = 2, 3, 4, 5, 6 . We also have 



T k (x)f (k \x) dx 

T k (x)dx = 0, , k = 2, . 

/ |T fc (x)| cfe = C fc , fc = 2 
Jo 



,5, 



,...,6, 



where 



C 2 

c fi 



19V19 

10125 : 
1 

604800 



C 3 



253 
360000 



C 4 



1 



14580 



1 



115200 



and 

where 
B 2 



1 

40 



max |T fc (x)| = B k 
xe[o,i] 



7 



fc = 2,3,4,5 



19V19 



20250 81000 



Ba 



5760 



58320 



(11) 



(12) 

(13) 
(14) 

(15) 
(16) 

(17) 
(18) 



Theorem 3 Let f G C fc (0, 1) and Zet r ) k , fre rea/ numbers such that 7^ < 
<T k , x e [0,1], fc = 2,3,4,5. Let S k = - /«(0) , fc = 

1,2,3,4. Then we have 



|iW)|<^^V 



/iy 2 ; . . . ? 5 



|i?fc(/)| < (^-i-7fe)5 fe , fc = 2,...,5 
< (T k -S k ^)B k , k = 2, 



,5, 



(19) 

(20) 
(21) 
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where Rk{f) o,re defined by (12), C k are defined by (15) and B k are defined 
by (18), fork = 2,3 ,4,5. 

Proof: Let C be an arbitrary constant. Then we have 

RkU)= I l T k (x)[fW(x)-C] dx= [ 1 T k (x)f^(x)dx, (22) 



for k = 2, 3, 4, 5 , since (13) holds. If we now choose C = (I\ + 7fc)/2 , then 
we get 



\Mf)\ 



TJx) 



dx 



< sup 

xe[o,i 



< 



f ik) (x) 

r fc - 7fc f 1 



Tfc + 7fe 



|Ti.(x)| dx 



T k (x)\ dx 



C k i 



(23) 



for k — 2, 3, 4, 5 . If we choose C = 7^ , then we have 



T fc (x) [/ (fc) (x)-7fc] dx 



< max |T fc (x)| / \f ik) {x)- lk \ dx 

= B k I [f {k \x)- lk ] dx 
Jo 

= (S k -i - lk)B k , 

for k = 2,3,4,5. 

In a similar way we can prove that (21) holds. 4k 

The estimations (19) are Peano-like bounds and they are generally (but 
not always) better than the usual Peano bounds. Namely, we know that T k (x), 
k — 2, ... ,6 are Peano kernels. The usual Peano error bounds are 



\Mf)\ < \\f 



Ml 



\T k (x)\ dx = C k \\f {k) \ 



fc = 2,...,6, 



where H/^H^ = sup x . e[01] |/ (fc) (x)| . If we choose j k = mi xe[0>1] f^(x) and 
T k = sup xg [ 01 ]/^ fc ^(x) , then Tk ~~ lk < Wf^W^- Thus, in this case, the error 
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bounds given by (19) are better than the Peano error bounds. In fact, they 
are equal if and only if r\ = — 7^ . This case (I\ = — 7^) is very rare in 
practice. 

Theoretically, we can derive better error bounds. Let us say something 
about the last assertion. We can verify that 

/ Tk(x)pj(x) dx — , j = 5 — k, k — 2, 3, 4, 
Jo 

where pj (x) is any polynomial of degree < j . Thus, 

f l T k (x) [fW{x)- Pj (x)] dx= I T k {x)f k \x)dx 
Jo Jo 

such that 

|^(/)|<||/ W -^LC fe , fc = 2,3,4, j = 5-k. 

The above estimations are theoretically better than the corresponding esti- 
mations in (19). 

We now give the above obtained results for an arbitrary interval [a,b]. 
The mapping x — (t — a)/(b — a) is a bijection from [a, b] onto [0, 1]. If we 
use this bijection then we find that the polynomials P k and Q k on the inter- 
val [a, b] have the forms: P Q (t) = 1 , Q (t) = 1 , A (*) = t - (23a + 76) /30 , 
Qi(t) = t — (7a + 236)/30 , etc. The polynomials can be also obtained by 
simple integration. For example, P2 can be obtained by integrating Pi and 
determining an additional constant such that (13) holds. (They are addi- 
tionally normalized.) We define the functions 

lk{t) - \ (-l) k Q k (t), te[^,b] ; [M) 

for k — 2, . . . , 6 . Here we choose a = , b = and use the notations 
introduced in Section 2. Using these functions we derive the following results. 

Corollary 4 Let f e C k (xj-i, Xj + \) , k = 2, 3, . . . , 6 . Then we have 

f(x)dx = ^ J —h-^ — h +R k {f), (25) 

Xj-l 

where 

R k (f)= n +1 f k (x)fW(x)dx,, (26) 

Jxj-! 

for k = 2, . . . , 6 and T k (x) defined by (24)- 
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Corollary 5 Let f G C k (xj_i, Xj + i) and let^ k , T k be real numbers such that 
Ik < f {k) {x) < T k , x G [xj-i, x j+1 ) , k = 2, 3, 4, 5 . Lei 



5Y 



k = 1,2,3,4 



JTien we /lave 



£*(/) <(S k - 1 -j k )E k h k+1 : 
R k (f) < (r k -S k ^)E k h k+1 



(27) 
(28) 
(29) 



where Rk(f) o,re defined by (24), D k = 2 k+1 C k (C k are defined by (15)) and 
E k = 2 k+1 B k (B k are defined by (18)), for /c = 2,3,4,5. 



Corollary 6 Let / G C 6 (xj_i, Xj+i) . Then we have 



Mf) 



<A>||/ (6) || h' 



(30) 



w/iere D 6 = 2 7 C 6 fC 6 defined by (16)) and WqW^ = sup^ j \g(x)\ 



4 Applications in numerical integration 

We define the partition tt = {a = xo < X\ < ■ ■ ■ < x 2n = b} of the inter- 
val [a, b] such that x i+ i = Xi + h , i — 0, 1, . . . , 2n — 1 , h = {b — a)/ (2n) . We 
also define the functions 

f, ( n / (-l) fe Pfci(:r) , x G ; 
kiW 1 (-l) fe g fc ,(x), X6(x 1)%1 ] ; 

for = 2, . . . , 6 , i = 1, 3, . . . , 2n — 1 , which correspond to the functions 
defined by (24) (for a = b = x i+ i). 

Theorem 7 Under the assumptions of Corollary 4 suppose that n and T k 
are given as above. Then we have 

b , 2n-l 

/ f(x)dx = - [7/i-i + 16/i + 7/ i+1 ] (31) 

10 odd 

[f'(b)-f(a)] + R k (f), 
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where 

£*(/)= E r +1 T ki {x)f^\x)dx, (32) 

j=l,j odd^j- 1 

for k = 2, . . . , 6 . 

Proof: Sum (25) over odd j from 1 to 2n — 1 to get (31-32), since 

2re-l 
j=l, j odd 

Theorem 8 Under the assumptions of Theorem 7 and Corollary 5 we have 

\R k (f)\<^p±D k h k (b-a), 
\R k {f)\<±{S k - 1 --y k )E k h k {b-a), 
\Mf)\ < \{V k -S k ^)E k h k {b-a), 

for k = 2, . . . , 5 . 

Proof: The proof follows immediately from Theorem 7 and Corollary 5. 

Theorem 9 Under the assumptions of Theorem 7 and Corollary 6 we have 

|i? 6 (/)|<^ 6 (fe-a)||/ (6) |L. 
Proof: The proof follows immediately from Theorem 7 and Corollary 6. 

Finally, let us compare the rule obtained in Theorem 7 with the standard 
composite Simpson's rule 

/b i 2n-l 

f(x)dx = - [/i-i + 4/i + /i+i] + JW) ■ ( 33 ) 

3=1,3 odd 

The terms 

h 2n_1 

3 E + 4/; + . 

3=1,3 odd 
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and 

, 2n-l 

15 E [7/i-i + Wi + 7/ i+ J , 

J=l, j odd 

require a same amount of calculations. The rule (31) has only one additional 
term with respect to the rule (33), namely 

lf'(b) - f'(a)} . 

Hence, the amount of calculations is approximately the same for both formu- 
lae. (Recall that function evaluations are generally considered the computa- 
tionally most expensive part of quadrature algorithms.) On the other hand, 
the rule (31) is exact for polynomials of degree < 5, while the rule (33) is 
exact for polynomials of degree < 3 . Furthermore, from Theorem 8 we have 

\Uf)\<^^h\b-a) : (34) 
while the standard estimation for the Simpson's rule is 

II f (4) ll 

\R s {J)\<^-^h\b-a). (35) 

Since, (34) is better than (35), the rule (31) has better approximation prop- 
erties than Simpson's rule. Thus, we expect that it will give better results in 
practice (in most cases). 



Example: here we show errors in estimating 

/ = Y^erf(l) = I exp(-x 2 )dx. (36) 

Using just two subintervals, h = 1/2, our formula (8) computes / ~ 0.746795 , 
whereas with four subintervals, that is h = 1/4, (8) gives I « 0.746824 which 
is correct to six decimal places. Figure 1 shows our rule converges quickly 
with decreasing grid spacing h, and is essentially exact to double precision 
with just 64 subintervals. 
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Figure 1: log-log plot of the errors of our integration rule (8), +'s, as a 
function showing the 0(/i 6 ) rate of convergence to the erf integral (36), 
compared to the 0(/i 4 ) convergence of the normal Simpson's rule, x's. 
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